function X_filtered = ecg_notch(X,F=50,Width=0.5,Depth=10,Samp_freq=550)
%% Notch filter for ECG
% Usage:
% X_filtered = ecg_notch(X);
% X_filtered = ecg_notch(X, F);
% X_filtered = ecg_notch(X, F, Width);
% X_filtered = ecg_notch(X, F, Width, Depth);
% X_filtered = ecg_notch(X, F, Width, Depth, Samp_freq);
%
% Inputs:
% X: Discrete time signal
% F: Notch in Hz
% Width: Width of the notch in Hz
% Depth: Depth of the notch in dB
% 
% Outputs:
% X_filtered: filtered signal
%
% Warnings:
% Do not overdesign the notch filter; filtfilt() may fail if the filter order is larger than 5
%
% Debug mode:
% If output variable is not specified, debug mode will automatically turned on.
% Filtered signal and Frequencncy response of the notch filter will be drawn in new figure window.
%
% Copyright by Group 8

Ws = [F-Width/2 F+Width/2] / (Samp_freq/2); % Stop band
Wp = [Ws(1)*0.7 Ws(2)*1.3]; % Pass band
Rp=1;  %Maximium ripple in the pass band (notch) in dB
Rs = Depth; %Minimium attenuation in the stop band (notch) in dB

[n,Wc] = buttord(Wp,Ws,Rp,Rs);

if(n>=5)
	warning('Filter order >=5, seems to overdesign bandreject filter!')
end

[b,a] = butter(n,Wc,'stop');

X_filtered = filtfilt(b,a,X);

%Debug mode:
if(!nargout)
figure(1)
freqz(b,a,512,Samp_freq);
figure(2)
t = [0:length(X)-1]/Samp_freq;
subplot(2,1,1), plot(t,X,';Raw signal;'),xlim([15 25])
ylabel('Voltage / V')
subplot(2,1,2), plot(t,X_filtered,';Filtered signal;'),xlim([15 25])
xlabel('time / s');
end

end


